﻿#####################################################
#													#
# Calculate P, etc from data file					#
# This script is an old branch of CalcPressure.py,	# 
# for generating the images in my thesis paper.		#
#													#
# Matt Ball (ball@brown.edu)						#
# 2009												#
#													#
#####################################################

import sys, Gnuplot, os, time
from time import strftime, localtime
from PyTrilinos import Epetra
from numpy import arange, std, array, zeros, ones, average
from scipy.integrate import quad, trapz
from math import pi, tanh

# for saving plots of (sortable) "scatterplot" functions, like amplitude & pressure rms
def saveArray(A, path, title, sortType, x_labels):
	g = Gnuplot.Gnuplot()
	# since enhanced mode interprets underscores as a control character
	enhancedTitle = title.replace("_","\\\_")
	#g('set title \"' + enhancedTitle + '\"')
	g('set term png size 1600,1200 enhanced font \"cmr12,24\"')
	g.xlabel(sortType)
	g.ylabel(title)
	#if not (sortType == None):
	#	title = title + "_by_" + sortType
	output = '\"~/Desktop/' + title + '.png\"'
	#output = output.replace("files","pressure")
	g('set output ' + output)
	
	if(x_labels != None):
		xtics_string = ""
		for i in range(len(x_labels)):
			xtics_string = xtics_string + '\"' + x_labels[i] + '\"' + str(i) + ', '
			if (i == len(x_labels)-1):
				# strip off the extra ', '
				xtics_string = xtics_string[:-2]
		
		g('set xtics (' + xtics_string + ')')
		g('set xtics nomirror out rotate')
	
	g.plot(A)
	
def plotStreamlines(A):
	title = "streamlines"
	g = Gnuplot.Gnuplot()
	
	numLines = 20 # the number of streamlines shown
	
	rows = A.M()
	cols = A.N()
	x = arange(-pi, pi, (2*pi)/rows)
	y = arange(0, p0, -abs(p0)/cols)
	g('set xrange [-pi:pi]')
	# since enhanced mode interprets underscores as a control character
	enhancedFileName = fileName.replace("_","\\\_")
	#g('set title \"' + title + ' (' + enhancedFileName + ')\"')
	#g.xlabel('q')
	#g.ylabel('height')
	g('set term png size 1600,1200 enhanced font \"cmr12,24\"')
	output = '\"~/Desktop/' + title + '_streamlines.png\"'
	#output = output.replace("files","pressure")
	g('set output ' + output)
	g('set data style dots')
	
	sls = zeros((rows*numLines,2))
	
	col = 0
	for i in range(rows*numLines):
		sls[i][0] = x[i%rows]
		sls[i][1] = A[i%rows][col]
		if (i % rows == 0):
			col += (cols/numLines)
	
	g.plot(sls)
	
def plotMat(A,title):
	# 'title' is "height", "pressure", etc.
	g = Gnuplot.Gnuplot()
	rows = A.M()
	cols = A.N()
	x = arange(-pi, pi, (2*pi)/rows)
	y = arange(0, p0, -abs(p0)/cols)
	g('set xrange [-pi:pi]')
	# since enhanced mode interprets underscores as a control character
	enhancedFileName = fileName.replace("_","\\\_")
	#g('set title \"' + title + ' (' + enhancedFileName + ')\"')
	g('set pm3d map')
	g.xlabel('q')
	g.ylabel('p')
	
	g('set term png size 1600,1200 enhanced font \"cmr12,24\"')
	
	output = '\"~/Desktop/' + title + '_2D.png\"'
	#output = output.replace("files","pressure")
	g('set output ' + output)
	
	g.splot(Gnuplot.GridData(A,x,y, binary=0))
	
def plotMatForAnimation(A, title, time):
	pass
	
def plotMat3D(A,title):
	# 'title' is "height", "pressure", etc.
	rows = A.M()
	cols = A.N()
	g = Gnuplot.Gnuplot()
	x = arange(-pi, pi, (2*pi)/rows)
	y = arange(0, p0, -abs(p0)/cols)
	g('set xrange [-pi:pi]')
	# since enhanced mode interprets underscores as a control character
	enhancedFileName = fileName.replace("_","\\\_")
	#g('set title \"' + title + ' (' + enhancedFileName + ')\"')
	g('set pm3d')
	g.xlabel('q')
	g.ylabel('p')
	g('unset surface')
	g('set view ,225,1,1')
	g('set ytics offset -.5,-.5')
	g('set term png size 1600,1200 enhanced font \"cmr12,24\"')
	
	output = '\"~/Desktop/' + title + '_3D.png\"'
	#output = output.replace("files","pressure")
	g('set output ' + output)
	
	g.splot(Gnuplot.GridData(A,x,y, binary=0))
	
# Trim edge entries for nicer printing
def trim(A):
	rows = A.M()
	cols = A.N()
	A_ = Epetra.SerialDenseMatrix(rows-2,cols-2,False)
	for i in range(rows-2):
		for j in range(cols-2):
			A_[i][j] = A[i+1][j+1]

	return A_

	
# make a directory called "pressure" whenever there is a directory called "files"
# ignore "pressure" if it's already there
def makePressureDirs(path):
	for dirpath, dirnames, filenames in os.walk(path):
		for dir in dirnames:
			if (dir == "files"):
				path = os.path.join(dirpath,"pressure")
				try:
					os.mkdir(path)
				except (OSError):
					pass
	
def readFile(file, name, saveIt):
	global p0
	global fileName
	global saveToDisk
	global Q
	global gamma
	global type
	global A # used for tanh
	global b
	
	fileName = name
	saveToDisk = saveIt
	
	type = file.readline() # exp or step
	rows = int(file.readline())
	cols = int(file.readline())
	p0 = float(file.readline())
	Q = float(file.readline())
	
	# the size of gamma here is *rows* since the matrix gets transposed after being read in
	if (type == "exp\n"):
		A = float(file.readline())
		b = float(file.readline())
		if (b == 0):
			type = "exp"
			gamma = A*ones(rows) # only for b=0 since gamma = A * exp[b*t] in general (for 'exp')
		else:
			# VERY SPECIAL CASE - tanh - NOT EXP at all!
			type = "tanh"
			gamma = zeros(rows) # to get delta_p in myGamma(...)
			
	elif(type == "step\n"):
		p1 = float(file.readline())
		p2 = float(file.readline())
		a0 = float(file.readline())
		a1 = float(file.readline())
		a2 = float(file.readline())
		type = "step"
		gamma = step_gamma(rows, p0, p1, p2, a0, a1, a2)
	
	file.readline() # blank line

	h = Epetra.SerialDenseMatrix(cols,rows,False)
	for i in range(rows):
		for j in range(cols):
			try:
				h[j][i] = float(file.readline())
			except ValueError:
				print "ValueError occurred somewhere near (probably after) line", `(i*cols)+j`
				print "The file is probably shorter than expected, sorry :("
				return
	
	h = centerWave(h)

	# height
	#plotMat3D(h, "height")
	#plotMat(h, "height")
	#plotMatForAnimation(h, "height", os.stat(fileName).st_mtime)
	
	# pressure
	pressure = P(h)
	plotMat3D(pressure, "pressure")
	#plotMat(pressure, "pressure")
	
	# streamlines
	#plotStreamlines(h)
	
	# P_x and P_y
	p_x = P_x(h, pressure)
	p_y = P_y(h, pressure)
	plotMat3D(trim(trim(p_x)),"P_x")
	plotMat3D(trim(trim(p_y)),"P_y")
	#plotMat(trim(trim(p_x)),"P_x")
	#plotMat(trim(trim(p_y)),"P_y")
	
# build a 1D array of \Gamma(p) for step - should be [a2, a2, a2, ..., a1, a1, a1, ..., a0, a0, a0]
def step_gamma(numElts, p0, p1, p2, a0, a1, a2):
	delta_p = abs(p0)/numElts
	temp = zeros(numElts)

	for i in range(int(abs(p2)/delta_p)):
		temp[i]= a2
	for i in range(int(abs(p2)/delta_p), int(abs(p1)/delta_p)):
		temp[i] = a1
	for i in range(int(abs(p1)/delta_p), int(abs(p0)/delta_p)):
		temp[i]= a0
	
	return temp

def myGamma(numElts,delta_p):
	Gamma = zeros(numElts)
	
	if (type == "exp") or (type == "step"):
		for i in range(numElts):
			Gamma[i] = trapz(gamma[:i+1], None, delta_p)
	
	elif (type == "tanh"):
		
		for i in range(numElts):
			# ExpWaveProb.h - line 29 for the formula (their Gamma is my ɣ)
			# 		// double Gamma(double p) { return A*tanh(B*(-abs(p)+1)); }
			Gamma[i] = quad(lambda p: A*tanh(b*(-abs(p)+1)), 0, delta_p*(i+1))[0]
	
	return Gamma
	
def centerWave(h):
	# the row at p = 0, iterate over q
	rows = h.M()
	cols = h.N()
	
	indexOfCenter = 0
	max = 0
	for i in range(rows):
		if (h[i][0] > max):
			max = h[i][0]
			indexOfCenter = i
	
	shiftBy = (rows/2) - indexOfCenter # this must be an int
	
	newWave = Epetra.SerialDenseMatrix(rows,cols,False)
	
	for i in range(rows):
		for j in range(cols):
			newWave[i][j] = h[(i-shiftBy)%rows][j]
	
	return newWave

# uses periodic boundary condition to evaluate centered difference at edges
def dq(h):
	rows = h.M()
	cols = h.N()
	
	delta_q = 2*pi/rows
	
	dq = Epetra.SerialDenseMatrix(rows,cols,False)
	for i in range(0,rows):
		for j in range(0,cols):
			# don't need to worry about negative indices, python takes care of them
			dq[i][j] = h[(i+1)%rows][j] - h[i-1][j]
			
	return dq/(2*delta_q)

# uses forward or backward difference at edges
def dp(h):
	rows = h.M()
	cols = h.N()
	
	delta_p = abs(p0)/cols
	
	dp = Epetra.SerialDenseMatrix(rows,cols,False)
	for i in range(0,rows):
		for j in range(0,cols):
			if (j==0):
				# forward
				dp[i][j] = (h[i][j+1] - h[i][j])/delta_p
				
			elif (j==cols-1):
				# backward
				dp[i][j] = (h[i][j] - h[i][j-1])/delta_p
				
			else:
				# general case
				dp[i][j] = (h[i][j+1] - h[i][j-1])/(2*delta_p)
			
	return dp
	
def P(h):
	rows = h.M()
	cols = h.N()
	
	hq = dq(h)
	hp = dp(h)
	
	g = 9.806 # For comparison to legacy code ONLY
	#g = 9.816
	
	delta_p = abs(p0)/cols
	Gamma = myGamma(cols, delta_p)
	
	P = Epetra.SerialDenseMatrix(rows,cols,False)
	for i in range(0,rows):
		for j in range(0,cols):
			P[i][j] = - .5*(1 + hq[i][j]**2)/(hp[i][j]**2) - g*h[i][j] - Gamma[j] + .5*Q
	
	return P
	
def P_x(h,P):
	rows = h.M()
	cols = h.N()
	
	hq = dq(h)
	hp = dp(h)
	
	Pq = dq(P)
	Pp = dp(P)
	
	Px = Epetra.SerialDenseMatrix(rows,cols,False)
	for i in range(rows):
		for j in range(cols):
			Px[i][j] = Pq[i][j] - (Pp[i][j])*(hq[i][j])/(hp[i][j])
	
	return Px

def P_y(h,P):
	rows = h.M()
	cols = h.N()
	
	hp = dp(h)
	
	Pp = dp(P)
	
	Py = Epetra.SerialDenseMatrix(rows,cols,False)
	for i in range(rows):
		for j in range(cols):
			Py[i][j] = Pp[i][j]/hp[i][j]
	
	return Py
	
# methods borrowed from Python v2.6 (kind of hacked)
def relpath(path, start='.'):
	
	# Return the longest prefix of all list elements.
	def commonprefix(m):
		"Given a list of pathnames, returns the longest common leading component"
		if not m: return ''
		s1 = min(m)
		s2 = max(m)
		for i, c in enumerate(s1):
			if c != s2[i]:
				return s1[:i]
		return s1
	
	"""Return a relative version of a path"""

	if not path:
		raise ValueError("no path specified")

	sep = '/'
	pardir = '..'
	
	start_list = os.path.abspath(start).split(sep)
	path_list = os.path.abspath(path).split(sep)

	# Work out how much of the filepath is shared by start and path.
	i = len(commonprefix([start_list, path_list]))

	rel_list = [pardir] * (len(start_list)-i) + path_list[i:]
	if not rel_list:
		return curdir
	return os.path.join(*rel_list)


# returns the standard deviation of the pressure along p=0
# this is not efficient since it computes the pressure over the entire wave! (optimize later)
# percentAroundCenter in [0,100]
def getPressureStdev(file, percentAroundCenter):
	global p0
	global type
	global gamma
	global type
	global Q
	
	type = file.readline() # exp or step
	rows = int(file.readline())
	cols = int(file.readline())
	p0 = float(file.readline())
	Q = float(file.readline())
	
	if (type == "exp\n"):
		A = float(file.readline())
		b = float(file.readline())
		type = "exp"
		gamma = A*ones(cols) # only for b=0 since gamma = A * exp[b*t] in general (for 'exp')
	elif(type == "step\n"):
		# confirm which is which
		p1 = float(file.readline())
		p2 = float(file.readline())
		a0 = float(file.readline())
		a1 = float(file.readline())
		a2 = float(file.readline())
		type = "step"
		gamma = step_gamma(cols, p0, p1, p2, a0, a1, a2)
	
	file.readline() # blank line
	
	temp = Epetra.SerialDenseMatrix(rows,cols,False)
	for i in range(rows):
		for j in range(cols):
			temp[i][j] = float(file.readline())

	# h
	h = Epetra.SerialDenseMatrix(cols,rows,False)
	for i in range(cols):
		for j in range(rows):
			h[i][j] = temp[j][i]

	h = centerWave(h)
	pressure = P(h)

	lowerBound = int(rows/2 - (rows*percentAroundCenter)/200)
	upperBound = int(rows/2 + (rows*percentAroundCenter)/200)
	numElts = upperBound - lowerBound

	topRow = zeros(numElts)
	for i in range(numElts):
		topRow[i] = pressure[i+lowerBound][0]
	
	return std(topRow)
	
def getAmplitude(file):
	type = file.readline() # exp or step
	rows = int(file.readline())
	cols = int(file.readline())
	p0 = float(file.readline())
	Q = float(file.readline())
	
	if (type == "exp\n"):
		A = float(file.readline())
		b = float(file.readline())
	elif(type == "step\n"):
		# p1, p2, a1, a2, a3
		file.readline()
		file.readline()
		file.readline()
		file.readline()
		file.readline()
	
	file.readline() # blank line
	
	# cols, not rows, because of the transpose issue
	topRow = zeros(cols)
	for i in range(cols):
		topRow[i] = file.readline()
	
	return max(topRow) - average(topRow)
	
def getQ(file):
	type = file.readline() # exp or step
	rows = int(file.readline())
	cols = int(file.readline())
	p0 = float(file.readline())
	Q = float(file.readline())
	
	return Q

def byDate():
	pressureRmsAndAmplitude("date")
	
def byQ():
	pressureRmsAndAmplitude("Q")

# amplitude and pressure RMS as functions of sortType (either date or Q)
def pressureRmsAndAmplitude(sortType):
	
	if (sortType != "date") and (sortType != "Q"):
		# in case the specified sort type isn't one of the implemented ones
		# check once at the beginning so I don't have to raise an exception later
		raise Exception, "sort type '" + sortType + "' not available"
	
	if ((sys.argv[2][-6:] == "files/") or (sys.argv[2][-5:] == "files")):
		print "Input directory must be a directory containing 'files' and not 'files' itself."
	else:
		makePressureDirs(sys.argv[2])
		numPoints = 0
	
		# count the number of data points we'll be generating
		for dirpath, dirnames, filenames in os.walk(sys.argv[2]):
			for f in filenames:
				if (dirpath[-6:] == "/files"):
					if (f[:4] == "sol."):
						numPoints = numPoints + 1
		
		print "Computing pressure RMS and amplitude values for", numPoints, "files (sorted by", sortType + ")"
		rms = zeros((numPoints,2))
		amp = zeros((numPoints,2))
		
		x_array = zeros(numPoints)
		x_labels = None
		
		if (sortType == "date"):
			x_labels = zeros(numPoints, dtype='|S20') # 20 = (max) length of the strings in the array
		
		i = 0
		
		# now actually do it
		for dirpath, dirnames, filenames in os.walk(sys.argv[2]):
			
			if (sortType == "date"):
				# this only works for date modified, not for date created
				filenames.sort(key = lambda f:os.stat(os.path.join(dirpath,f)).st_mtime)
			#else (sortType == "Q"):
				# filenames is already in alphabetical order, so do nothing
			#	pass
			
			for f in filenames:
				
				if (dirpath[-6:] == "/files"):
					if (f[:4] == "sol."):
						name = os.path.join(dirpath, f)
						
						print "(" + str(i+1) + ")\t" + name
						
					#	file = open(name,"r")
					#	
					#	if(len(sys.argv) == 4):
					#		rms[i][1] = getPressureStdev(file, int(sys.argv[3]))
					#	else:
					#		rms[i][1] = getPressureStdev(file, 100)
					#	
					#	file.close()
						file = open(name,"r")
						amp[i][1] = getAmplitude(file)
						file.close()
						
						if(sortType == "Q"):
							file = open(name,"r")
							x_array[i] = getQ(file)
							file.close()
						else:
							# sortType equals "date"
							x_array[i] = i
							x_labels[i] = str(strftime("%Y-%m-%d %H:%M:%S",localtime(os.stat(os.path.join(dirpath,f)).st_mtime)))
							
						i += 1
						
		for j in range(numPoints):
		#	rms[j][0] = x_array[j]
			amp[j][0] = x_array[j]
		
		#pressureTitle = ""
		#if (len(sys.argv) == 4):
		#	pressureTitle = "pressure_rms_" + sys.argv[3] + "%"
		#else:
		#	pressureTitle = "pressure_rms_100%"
		#
		#saveArray(rms, dirpath, pressureTitle, sortType, x_labels)
		saveArray(amp, dirpath, "amplitude", sortType, x_labels)

# creates/updates the XML database (<waveset>...</waveset>) with entries of the form:
#	<wave class="exp" p0="-2" Q="21.5" amplitude="3.4" path="foo/bar/baz" A="5"
#		resolution_q="500" resolution_p="250" >
#		<graph dims="2" type="pressure" path="foo/bar/graph" />
#	</wave>
# important: path names are relative to the directory this script was called from
# only operates on 'sol.' files found in directories named 'files'
def makeDatabase():
	
	print "\nMake sure to run this script from the wave data root directory when generating the database.\n"
	
	dbName = "database.xml"
	dbPath = os.path.join(sys.argv[2], dbName)
	
	file = None
	
	if (os.path.exists(dbPath)):
		# there already exists a file called database.xml in the directory spec'd by argv[2]
		YorN = raw_input("Overwrite old database file? ")
		yorn = YorN.lower()
		if (yorn not in ["y", "yes"]):
			print "exiting..."
			return
		else:
			pass
			
	# check that argv[2] is a directory
	if (os.path.isdir(sys.argv[2])):
		# it is
		file = open(dbPath, "w") # w mode overwrites the old file, w+ appends
		
		file.write("<waveset>\n")
		numEntries = 0
		numFailedEntries = 0
		
		# Write the entries
		for dirpath, dirnames, filenames in os.walk(sys.argv[2]):

			for f in filenames:
						
				if (dirpath[-6:] == "/files"):
				
					if (f[:4] == "sol."):
						name = os.path.join(dirpath, f)
						
						thisEntry = dbEntry(name)
						# don't write empty entries
						if not (thisEntry == "\t<wave \n\t</wave>\n"):
							file.write(thisEntry)
							numEntries += 1
						else:
							numFailedEntries += 1
			
		file.write("</waveset>")
		file.close()
		print "database.xml is complete, with", numEntries, "waves successfully entered and", numFailedEntries, "waves failed"
		print "Currently, this file must be broken up into db_exp.xml, db_step.xml, and db_tanh.xml for use by the web app."
		return
		
	elif (os.path.isfile(sys.argv[2])):
		# it's a file, so we'll treat it as a single wave to operate on
		print dbEntry(sys.argv[2])
		return
	else:
		print "Invalid path, exiting..."
		return

# returns a string representing a database entry for a single file
def dbEntry(filePath):

	file = open(filePath,"r")
	amp = getAmplitude(file)
	file.close()
	
	entry = "\t<wave "

	try:
		file = open(filePath)
		time = int(os.stat(filePath).st_mtime)
		type = file.readline() # exp or step
		rows = int(file.readline())
		cols = int(file.readline())
		p0 = float(file.readline())
		Q = float(file.readline())
		
		# the size of gamma here is *rows* since the matrix gets transposed after being read in
		if (type == "exp\n"):
			A = float(file.readline())
			b = float(file.readline())
			if (b == 0):
				type = "exp"
				entry += "class=\"" + type + "\" p0=\"" + str(p0) + "\" Q=\"" + str(Q)
				entry += "\" amplitude=\"" + str(amp) + "\" path=\"" + relpath(filePath)
				entry += "\" A=\"" + str(A) + "\" resolution_q=\"" + str(rows)
				entry += "\" resolution_p=\"" + str(cols) + "\" time=\"" + str(time) + "\">"
				
			else:
				# VERY SPECIAL CASE - tanh - NOT EXP at all!
				type = "tanh"
				entry += "class=\"" + type + "\" p0=\"" + str(p0) + "\" Q=\"" + str(Q)
				entry += "\" amplitude=\"" + str(amp) + "\" path=\"" + relpath(filePath)
				entry += "\" A=\"" + str(A) + "\" B=\"" + str(b) + "\" resolution_q=\"" + str(rows)
				entry += "\" resolution_p=\"" + str(cols) + "\" time=\"" + str(time) + "\">"
				
		elif(type == "step\n"):
			p1 = float(file.readline())
			p2 = float(file.readline())
			a0 = float(file.readline())
			a1 = float(file.readline())
			a2 = float(file.readline())
			type = "step"
			entry += "class=\"" + type + "\" Q=\"" + str(Q) + "\" p0=\"" + str(p0)
			entry += "\" p1=\"" + str(p1) + "\" p2=\"" + str(p2) + "\" a0=\"" + str(a0)
			entry += "\" a1=\"" + str(a1) + "\" a2=\"" + str(a2) + "\" amplitude=\"" + str(amp)
			entry += "\" path=\"" + relpath(filePath) + "\" resolution_q=\"" + str(rows)
			entry += "\" resolution_p=\"" + str(cols) + "\" time=\"" + str(time) + "\">"
		
		file.close()
	
	except:
		pass
	
	entry += "\n\t</wave>\n"
	return entry

# only operates on 'sol.' files found in a directory named 'files'
def batch():
	# count the number of data points we'll be generating
	numPoints = 0
	for dirpath, dirnames, filenames in os.walk(sys.argv[2]):
		for f in filenames:
			if (dirpath[-6:] == "/files"):
				if (f[:4] == "sol."):
					numPoints = numPoints + 1
		
	print "Generating plots for", numPoints, "waves. This could take a while."
	YorN = raw_input("Proceed? ")
	yorn = YorN.lower()
	if (yorn not in ["y", "yes"]):
		print "exiting..."
		return
	else:
		pass

	makePressureDirs(sys.argv[2])
	
	# ugly hack to periodically print out % complete. add more elements to this array for more frequent messages.
	hacked = array([True, True, True, True, True, True, True, True, True, True])
	percentComplete = 0
	
	for dirpath, dirnames, filenames in os.walk(sys.argv[2]):
	
		for f in filenames:
					
			if (dirpath[-6:] == "/files"):
			
				if (f[:4] == "sol."):
					try:
						percentComplete += (1./numPoints)*100
						
						for i in range(len(hacked)):
							if (percentComplete >= (i*(100/len(hacked)))) and hacked[i]:
								print str(int(percentComplete))+"% complete"
								hacked[i] = False
								
						name = os.path.join(dirpath, f)
						print name
						file = open(name,"r")
						readFile(file, name, True)
						file.close()
					except:
						pass
							
def single():
	name = sys.argv[2]
	file = open(name,"r")
	readFile(file, name, False)
	file.close()
	
def main():

	funcmap = dict(single=single, batch=batch, pRMS_amp_byQ=byQ, pRMS_amp_byDate=byDate, db=makeDatabase)
	def usage():
		print 'Usage: python %s [%s]' % (sys.argv[0], ' | '.join(funcmap.keys())) + ' path [RMS % around center]'
	funckey = (len(sys.argv) >= 3) and (len(sys.argv) <= 4) and (sys.argv[1]) or (None)
	func = funckey in funcmap and funcmap[funckey] or usage	
	func()
	

if __name__ == '__main__':
    main()